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Abstract 

The  solution  of  nonlinear  optimal  control  problems  and  the  computation 
of  optimal  control  laws  is  a  difficult  task  that  requires  a  fast  response 
in  many  real  applications.  During  the  three  years  of  the  AFOSR  grant 
we  attacked  the  problem  via  the  numerical  solution  of  Hamilton-Jacobi 
equations  describing  the  value  function.  In  the  framework  of  viscosity  so¬ 
lutions,  the  value  function  is  uniquely  characterized  and  it  allows  for  the 
computation  of  optimal  controls  in  feedback  form.  Unfortunately,  from  a 
computational  point  of  view,  solving  real  problems  in  high  dimension  via 
Hamilton-Jacobi  equations  is  still  a  huge  task  and  efficient  algorithms  are 
required.  Our  research  has  been  devoted  to  the  development  of  new  Fast 
Marching  methods  and  domain  decomposition  techniques  (or  a  combina¬ 
tion  of  them)  in  order  to  improve  CPU  times,  avoid  useless  computations 
and  reduce  memory  requirements. 


1  Introduction 

The  computation  of  optimal  controls  is  one  of  the  primary  goals  in  control  the¬ 
ory.  The  classical  approach  based  on  Pontryagin’s  Maximum  Principle  leads  to 
a  two  point  boundary  value  problem  for  a  system  of  ordinary  differential  equa¬ 
tions  describing  the  dynamics  and  the  co-state  and  provide  necessary  condi¬ 
tions  for  optimality.  Although  some  numerical  methods  have  been  proposed  for 
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that  problem,  it  is  well  known  that  this  approach  is  limited  by  the  following 
difficulties:  a)  the  initialization  of  the  co-state  equations  to  start  the  numerical 
procedure  can  be  very  difficult  and  often  requires  weeks  to  be  solved,  b)  the 
optimal  control  which  is  obtained  is  open-loop  and  c)  the  optimal  state-costate 
couple  just  satisfies  necessary  conditions  so  there  is  no  guarantee  that  it  corre¬ 
sponds  to  a  global  minimum  for  the  cost  functional  on  the  space  of  admissible 
controls.  Despite  those  limitations,  the  approach  based  on  the  Pontryagin’s 
Maximum  Principle  has  been  up  to  now  the  only  approach  which  has  produced 
feasible  solutions  for  real  industrial  problems. 

The  second  classical  approach  to  the  solution  of  optimal  control  problems 
is  Dynamic  Programming,  introduced  by  R.  Bellman  in  the  60’s  [7].  In  this 
approach,  a  central  role  is  played  by  the  value  function  u  of  the  problem,  defined 
as 

u(x)  =  inf  J(x,a(-)) 

a(-)eA 

and  the  reconstruction  of  control  is  done  starting  from  the  Bellman  equation 
associated  to  the  problem  once  the  value  function  has  been  found  (here  x  will 
represent  the  initial  position  of  the  system,  A  will  be  the  set  of  admissible  con¬ 
trols  and  J  is  the  cost  functional). 

The  development  of  the  theory  of  viscosity  solutions  in  the  last  twenty  years 
has  finally  given  a  general  framework  for  the  characterization  of  the  value  func¬ 
tion  in  all  the  classical  problems  of  deterministic  and  stochastic  control  theory 
(finite  horizon,  infinite  horizon,  minimum  time,  optimal  switching,  impulsive 
control,  pursuit-evasion  games).  Indeed,  under  very  general  assumptions,  one 
is  able  to  prove  that  the  unique  viscosity  solution  of  the  Bellman  (resp.  Isaacs) 
equation  associated  to  the  control  (resp.  game)  problem  is  the  value  function. 
The  advantage  of  this  approach  is  twofold.  It  gives  a  precise  characterization  of 
the  global  optimum  for  the  problem  and  it  allows  to  obtain  optimal  controls  in 
feedback  form.  This  has  motivated  a  large  research  activity  on  the  development 
of  numerical  methods  to  solve  the  Bellman  (Isaacs)  equation  as  well  as  on  the 
numerical  synthesis  of  feedback  controls  based  on  the  knowledge  of  an  approxi¬ 
mate  value  function. 

The  main  drawback  of  this  approach  is  due  to  the  well  known  “curse  of  di¬ 
mensionality”  of  Dynamic  Programming.  Indeed,  the  numerical  solution  of  a 
nonlinear  partial  differential  equation  in  high-dimension  is  a  difficult  task  and 
requires  new  ideas  and  algorithms.  This  project  had  the  aim  of  developing  new 
methods  coupling  the  semi-Lagrangian  (or  Dynamic  Programming)  approxima¬ 
tion  with  Fast  Marching  and  domain  decomposition  techniques  in  order  to  solve 
optimal  control  and  pursuit-evasion  games. 

2  Background 

In  this  section  we  briefly  review,  for  the  reader’s  convenience,  some  classical 
topics  underlying  the  new  methods  we  developed  during  the  three  years  of  the 
AFOSR  grant. 
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2.1  Semi-Lagrangian  schemes  for  control  problems  and 
games  via  Dynamic  Programming  equations 

The  Dynamic  Programming  approach  is  based  on  the  characterization  of  the 
value  function  of  the  problem  via  the  Bellman  equation,  a  nonlinear  first  order 
partial  differential  equation  for  deterministic  optimal  control  and  game  prob¬ 
lems.  For  instance,  the  Hamilton- Jacobi-Bellman  (HJB)  equation  for  minimum 
time  problems  (with  target  set  T)  is 

f  max{-/(r,  a)  •  V'u(x)}  —  1  =  0  x  G  Mn\T, 

<  (1) 

|  u(x )  =  g(x)  x  G  dT . 

This  leads  to  the  characterization  of  the  optimal  control  in  feedback  form,  which 
requires  the  exact  (or  approximate)  knowledge  of  u : 

a* hr)  =  argmin{/(x,  a)  •  Wu(x)}  ,  x  G  Rn\T  -  (2) 

aeA 

Similarly,  the  Hamilton- Jacobi-Isaacs  (HJI)  equation  for  pursuit-evasion  games 
is 

J  minma x{— /(x,  a,  b )  •  V^(x)}  —  1  =  0  x  G  Mn\T 

|  u(x)  =  ^(x)  x  G  dT 

and  the  corresponding  optimal  feedback  strategies  can  be  computed  from  u  by 

(a*(x),  b*(x))  =  argmaxmin{/(x,  a,  6)  •  Tu(x)}  ,  x  G  Mn\T.  (4) 
beB  aeA 

Let  us  introduce  the  basic  concepts  of  Semi-Lagrangian  (SL)  schemes  that  we 
will  use  as  a  building  block  for  our  methods. 

Consider  a  structured  grid  G  and  denote  its  nodes  by  i  =  1, ...  ,7V.  The 
space  step  is  assumed  to  be  uniform  and  equal  to  Ax  >  0. 

The  HJB  equation  (1)  can  be  discretized  by  means  of  a  discrete  Dynamic  Pro¬ 
gramming  Principle,  so  that  the  relationship  with  the  optimal  control  framework 
is  never  lost.  Standard  arguments  [5]  lead  to  the  following  scheme  (similar  ar¬ 
guments  apply  for  the  HJI  equation  (3)): 


u(xi)  =  min  <  u(xijCL)  + 


1  ^ i,a)  T  \  |  (  5  W 

aeA  [  \j(xi,a)\  J 

where  is  a  non- grid  point,  obtained  by  integrating,  until  a  certain  final  time 
r,  the  ordinary  differential  equation 

f  y(t)  =  f(y,a ),  te[o,r]  f  ) 

I  y( o)  =  x i 


and  then  setting  =  y(r).  To  make  the  scheme  fully  discrete,  the  set  of 

admissible  controls  A  is  discretized  in  Nc  points  and  the  minimum  in  (5)  is 
computed  by  a  direct  comparison  of  the  Nc  values. 

Depending  on  the  choice  of  r,  on  how  equation  (6)  is  solved  and  on  how  u(x^a) 
is  interpolated,  we  get  different  versions  of  the  semi-Lagrangian  scheme  (5). 
Here  we  list  two  different  discretizations  that  often  appear  in  the  literature: 
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•  2-points  SL.  This  scheme  is  used,  for  example,  in  [26]  and  [27].  Equation 
(6)  is  solved  by  an  explicit  forward  Euler  scheme  until  the  solution  inter¬ 
cepts  the  line  connecting  two  neighbouring  points  and  X{ ,2  (see  Figure 
la).  The  value  u(xi)  is  computed  by  linear  interpolation  of  the  values 
u(xi: i)  and  u(xi^)- 

•  3-points  SL.  This  scheme  is  used,  for  example,  in  [11].  Equation  (6)  is 
solved  by  an  explicit  forward  Euler  scheme  until  the  solution  is  at  distance 
Ax  from  where  it  falls  inside  the  triangle  of  vertices  x^  i,  x^2,  and  x^ 3 
(see  Figure  lb).  The  value  u(xi)  is  computed  by  linear  interpolation  of 
the  values  u(x^  1),  u(x^ 2)  and  u(x^ 3). 


I  ^,3|  ^  hcj,2 
;  /(a?i,oy  '  \ 


t  -r  t 

(b) 


Figure  1:  (a)  2-points  SL  scheme,  (b)  3-points  SL  scheme. 


2.2  Fast  Marching  Method  for  the  Eikonal  equation 

The  Fast  Marching  Method  (FMM)  has  been  introduced  by  Sethian  [25]  in  1996 
as  a  fast  solver  for  the  eikonal  equation: 

(  c(x)\\7u(x)\  =  1  xeRn\T , 

\  u(x)  =  g(x)  x  e  <9T, 

where  c,  g  are  Lipschitz  functions  and  the  target  set  T  is  a  closed  subset. 
Unlike  more  classical  fixed-point  methods  that  simply  iterate  “brute  force”  on 
the  whole  grid  until  convergence,  the  main  idea  of  FMM  is  to  visit  grid  nodes 
in  a  dynamic-dependent  order  exploiting  the  hyperbolic  nature  of  the  equation, 
namely  trying  to  follow  the  characteristic  flow  emanating  from  the  target  set  T- 
To  this  end  grid  nodes  are  iteratively  divided  in  three  subsets:  Accepted  (ACC) 
region,  Narrow  Band  ( NB )  region  and  Far  (FAR)  region  (see  Figure  2).  Nodes 
in  ACC  are  definitively  computed,  nodes  in  NB  are  computed  but  their  values 
are  not  yet  final  and  nodes  in  FAR  are  not  yet  computed. 

At  the  beginning  ACC  only  consists  of  nodes  belonging  to  the  target  set  T, 
whose  values  are  imposed  using  the  Dirichlet  condition  given  by  g.  Then  the  al¬ 
gorithm  finds  itself  a  correct  order  for  processing  all  the  remaining  nodes.  That 
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Figure  2:  Partition  of  grid  nodes  in  Accepted ,  Narrow  Band  and  Far. 


order  is  found  in  such  a  way  a  causality  principle  is  respected,  i.e.  the  value  of 
a  node  is  labelled  as  “final”  only  if  it  depends  on  already  “final”  values  and  it 
cannot  be  affected  by  the  future  computation. 

In  practice,  the  node  in  NB  with  the  minimal  value  enters  ACC  and  its  first 
neighbours  are  computed  by  the  underlying  scheme  entering  NB  (if  not  yet  in). 
We  remark  that  this  minimal  value  rule  corresponds  to  compute  the  solution  of 
the  equation  following  its  level  sets.  Indeed,  suppose  for  simplicity  that  g  =  0, 
so  that  <9T  is  the  zero  level  set  of  u.  Then,  the  node  in  NB  with  minimal  value 
identifies  the  minimal  variation  of  the  solution,  hence  the  direction  of  its  gradi¬ 
ent.  It  follows  that  NB  expands  under  the  gradient  flow  of  the  solution  itself, 
which  is  exactly  equivalent  to  say  that  NB  is,  at  each  step,  an  approximation 
of  a  level  set.  In  the  case  of  the  Eikonal  equation,  the  gradient  of  the  solution 
coincides  with  the  characteristic  field  of  the  equation,  hence  FMM  computes 
the  correct  solution,  namely  the  same  solution  obtained  with  a  full  grid  method 
employing  the  same  scheme. 

More  in  general,  FMM  still  works  for  problems  with  mild  anisotropy,  where 
gradient  lines  and  characteristics  define  small  angles  and  lie,  at  each  point,  in 
the  same  simplex  of  the  underlying  grid. 

On  the  other  hand,  when  a  strong  anisotropy  comes  into  play,  as  for  the 
anisotropic  Eikonal  equation  studied  in  [26],  FMM  fails  and  there  is  no  way 
to  compute  the  solution  following  its  level  sets. 

FMM  is  an  example  of  what  is  called  a  local  single-pass  scheme,  where  local 
refers  to  the  fact  that  the  computation  of  a  grid  point  involves  only  the  values 
of  neighbouring  points,  the  region  NB  is  only  one-cell  or  two-cell  large  and  no 
information  coming  from  FAR  region  is  used.  On  the  other  hand  single-pass 
means  that  each  grid  point  is  re-computed  at  most  r  times,  where  r  only  de¬ 
pends  on  the  considered  equation  and  the  grid  structure,  not  on  the  number  of 
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grid  points. 


2.3  Domain  Decomposition  for  Hamilton- Jacobi  equations 

The  Domain  Decomposition  is  a  classical  and  quite  general  technique  to  solve 
a  problem  on  a  huge  grid  by  dividing  it  into  sub-problems  of  a  reasonable  size. 
The  idea  is  rather  simple  and  intuitive:  the  domain  of  computation  is  split  in 
several  sub-domains  and  each  sub-domain  is  assigned  to  a  processor,  so  that 
the  whole  computation  is  carried  on  in  parallel  (see  for  instance  [24]).  For  an 
application  of  this  technique  in  the  context  of  Hamilton- Jacobi  equations  we 
refer  the  reader  to  [12,  18].  The  delicate  point  for  the  efficiency  of  this  kind  of 
parallel  methods  is  the  fact  that  processors  must  share  some  parts  of  their  sub- 
domains  (at  least  the  sub-domain  boundaries)  and  must  communicate  during 
the  computation  using  transmission  conditions  (see  Figure  3). 


Q 

2 

■ 

r 

Q 

Q 

3 

4 

(b) 


Figure  3:  Classical  domain  decomposition  with  overlapping  boundaries:  (a)  2 
sub-domains,  (b)  4  sub-domains. 

Since,  in  general,  there  is  no  a-priori  knowledge  on  where  the  correct  information 
will  flow  from,  due  to  the  arbitrariness  of  the  decomposition,  these  communica¬ 
tions  must  take  place  at  every  iteration.  This  is  a  real  bottleneck  for  the  parallel 
computation,  since  each  processor  continues  working  until  a  global  convergence 
condition  is  satisfied. 


3  Objectives 

Let  us  recall  the  original  objectives  of  this  research  project: 

3.1  First  year 

Analysis  and  development  of  Fast  Marching  methods  for  first  order  Hamilton- 
Jacobi  equation  with  convex  Hamiltonians.  Applications  to  optimal  control 
problems  in  low  dimension. 
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3.2  Second  year 

Extension  of  Fast  Marching  methods  to  high-dimensional  control  problems. 
Analysis  of  convergence,  a-priori  error  estimates.  Implementation  of  the  al¬ 
gorithms  and  coupling  with  the  domain  decomposition  technique.  Parallel  al¬ 
gorithms. 

3.3  Third  year 

Extension  of  Fast  Marching  methods  to  non-convex  Hamiltonians  and  applica¬ 
tions  to  differential  games.  Analysis  of  convergence,  a-priori  error  estimates. 
Implementation  of  the  algorithms  and  coupling  with  the  domain  decomposition 
technique.  Parallel  algorithms. 

4  Findings 

The  results  that  have  been  obtained  during  this  project  are  mainly  related  to 
the  development  of  new  ideas  and  algorithms.  Only  in  few  cases  we  obtained  a 
complete  analysis  of  the  new  methods,  giving  new  bounds  for  their  complexity 
and  acceleration.  However,  it  is  interesting  to  note  that  the  general  bounds 
for  global  semi-lagrangian  methods  (see  e.g.  [13,  14])  still  apply.  The  extensive 
testing  of  our  methods  on  a  large  number  of  benchmarks  shows  that  the  methods 
are  really  effective  and  improve  the  results  in  the  literature. 

4.1  A  Parallel  Fast  Marching  Method  for  the  Eikonal  equa¬ 
tion 

The  classical  FMM  described  in  Section  2.2  works  in  a  highly  sequential  way, 
due  to  the  fact  that  the  minimal  value  rule  accepts  only  one  grid  node  in  NB 
at  a  time.  This  is  the  main  reason  explaining  the  scarce  literature  on  paral¬ 
lelization  techniques  for  FMM.  To  our  best  knowledge,  only  few  contributions 
are  available  on  this  topic  (see  for  instance  [20,  28,  29])  dealing  with  parallel 
algorithms  based  on  domain  decomposition  methods,  in  some  cases  optimized 
for  Graphics  Processing  Units. 

In  this  direction  we  developed  a  parallelization  approach  for  FMM  which  is  not 
based  on  a  domain  decomposition  technique  and  it  is  also  easier  to  implement 
than  Hermann’s  method  proposed  in  [20]. 

The  main  idea  is  to  split  among  threads  only  the  part  of  the  boundary  of  the 
domain  where  the  Dirichlet  conditions  are  given,  rather  than  splitting  the  whole 
computational  domain.  The  sub-domains  resulting  from  this  splitting  are  then 
used  by  each  thread  as  starting  points  for  computing  independent  solutions. 
Thereby,  the  interactions  among  threads  are  realized  by  a  relatively  simple  pro¬ 
cedure  which  guarantees  the  correct  order  of  acceptance  and  keeps  the  overlap¬ 
ping  areas  to  the  minimum.  In  the  end,  the  computational  domain  is  automat¬ 
ically  divided  among  threads  in  a  way  which  depends  on  the  dynamics  of  the 
problem.  Several  tests  showed  that  our  method  is  computationally  competitive 
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with  Hermann’s  method,  in  some  cases  even  more  efficient. 

Figure  4  shows  a  comparison  of  the  speed  up  for  our  algorithm  and  Hermann’s 
method  applied  to  the  Eikonal  equation. 


Figure  4:  Speedup  comparison. 

In  this  test  the  domain  is  a  square  and  boundary  conditions  are  imposed  on 
a  small  circle  located  in  a  corner  of  the  domain.  The  pre-determined  domain 
decomposition  employed  by  Hermann  splits  the  domain  into  equal  parts.  This 
choice  is  clearly  not  optimal,  since  at  the  beginning  of  the  computation  only  the 
thread  assigned  to  the  sub-domain  containing  the  circle  is  really  active.  Our 
algorithm  is  able  to  assign  a  job  to  each  core  from  the  very  beginning  and  never 
degenerates  into  a  serial  mode.  We  refer  the  reader  to  [PI]  for  further  details. 

4.2  Progressive  Fast  Marching  Method  (PFMM)  for  con¬ 
trol  problems  and  games 

Several  efforts  have  been  made  to  extend  the  classical  FMM,  in  order  to  increase 
its  efficiency  and  manage  more  general  equations/problems  (see  e.g.  [22]  and 
the  references  therein).  As  described  in  Section  2.2,  FMM  is  designed  on  the 
Eikonal  equation,  exploiting  the  fact  that  the  characteristic  field  of  the  equation 
coincides  with  the  gradient  of  the  solution  itself.  This  property  does  not  hold  for 
general  Hamilton- Jacobi  equations,  hence  FMM  fails  in  computing  the  correct 
solution,  especially  in  problems  with  strong  anisotropy.  To  our  opinion,  the  most 
valuable  work  aiming  at  overcoming  this  difficulties  is  contained  in  the  paper 
[26],  where  an  Ordered  Upwind  Method  has  been  proposed  to  solve  anisotropic 
equations  in  a  Fast  Marching  fashion.  This  method  produces  good  results  for 
several  difficult  tests,  but  the  spirit  of  FMM  is  somehow  betrayed,  in  the  sense 
that  the  locality  property  of  the  method  is  completely  lost.  Moreover,  the 
computed  solution  does  not  coincide  with  the  solution  obtained  by  a  full  grid 
method  using  the  same  scheme.  The  error  depends  on  the  anisotropy  of  the 
problem  and  vanishes  only  in  the  limit,  as  the  space  step  of  the  grid  goes  to 


zero. 


In  this  direction  we  investigated  the  possibility  of  keeping  the  locality  and 
single-pass  properties  of  FMM,  but  replacing  the  minimal  value  rule  as  accep¬ 
tance  criterium  for  nodes  in  NB. 

In  the  framework  of  optimal  control,  we  analyzed  the  behavior  of  the  approxi¬ 
mate  value  function  solving  the  HJB  equation  (1)  subject  to  local  perturbations, 
aiming  at  understanding  how  NB  nodes  depend  on  the  already  accepted  ones. 
We  observed  that,  over  a  certain  threshold  of  the  perturbation,  only  a  small 
number  of  grid  points  of  the  narrow  band  do  not  vary,  and  this  suggests  that 
they  have  to  be  processed  before  the  others. 

To  be  more  precise,  consider  the  situation  sketched  in  Figure  5a.  Red  nodes  are 
in  ACC  ,  green  nodes  are  in  NB  and  black  nodes  are  in  FAR. 


Figure  5:  Dependency  of  NB  nodes  on  local  perturbations  of  the  value  function. 

For  each  x  G  NB  consider  its  first  neighbors,  namely  the  nodes  needed  to  com¬ 
pute  the  value  u(x).  The  arrow  represents  the  optimal  dynamics  at  x  in  the 
HJB  equation  (1).  Now  replace  the  values  of  all  the  not- accepted  neighbors  of 
x  by  a  tentative  value  u  (blue  nodes  in  Figure  5b).  A  small  tentative  value  can 
change  the  optimal  control  and  then  the  value  u(x)  (see  Figure  5c).  But,  if  we 
progressively  increase  this  tentative  value,  it  can  happen  that,  over  a  certain 
threshold  (see  light  blue  nodes  in  Figure  5d),  the  “new”  optimal  control  is  no 
longer  favorable  and  the  scheme  comes  back  to  choose  the  “original”  optimal 
control  which  only  depends  on  accepted  nodes.  In  this  way,  we  actually  mimic 
the  influence  that  still  unknown  values  (NB  and  FAR  neighbors)  can  have  on 
u(x)  in  the  future.  Then,  we  define  the  threshold  a(x)  as  the  minimal  tentative 
value  such  that  u(x)  does  not  change  for  any  u  >  a(x).  This  threshold  quanti¬ 
fies  the  dependency  of  u(x)  on  its  not-accepted  neighbors  and  then  it  appears 
natural  to  consider  the  node  in  NB  with  minimal  threshold  as  the  most  correct. 
This  argument  produces  a  new  criterium  for  ordering  the  grid  nodes  (different 
from  the  minimal  value  rule),  showing  the  ability  to  compute  the  correct  solu¬ 
tion  for  a  wide  class  of  Hamilton- Jacobi  equations,  also  in  the  case  of  nonconvex 
Hamiltonians,  e.g.  for  Hamilton- Jacobi-Isaacs  equations  related  to  Pursuit- 
Evasion  games.  Moreover,  the  Progressive  Fast  Marching  Method  (PFMM)  is  a 
coherent  extension  of  the  Fast  Marching  method,  since  it  gives  back  the  classical 
FMM  when  applied  to  isotropic  equations  of  eikonal  type.  Indeed,  in  these  cases 
the  node  with  minimal  threshold  is  also  the  node  with  minimal  value.  We  refer 
the  reader  to  [P2]  for  further  details. 

Figure  6  presents  the  result  of  PFFM  compared  with  the  classical  FMM  in  solv- 
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ing  the  anisotropic  Eikonal  equation,  that  can  be  re-casted  as  an  HJB  equation 
by  choosing  in  (1) 


f(x,y,a)  =  a(l  +  (\a1+ida2)2)  2 

H,  A  >  0,  a  =  (oi,a2)  e  A  =  B(0, 1),  T  ={(0,0)}. 

In  Figure  6a-b  we  show  the  level  sets  of  the  wrong  solution  computed  by  FMM 
against  the  correct  solution  computed  by  PFFM,  in  Figure  6c-d  we  show  the 
corresponding  orders  of  acceptance  at  an  intermediate  stage. 


Figure  6:  Anisotropic  Eikonal  equation.  Level  set  of  the  solution:  (a)  FMM, 
(b)  PFFM.  Order  of  acceptance:  (c)  FMM,  (d)  PFFM. 


Figure  7  presents  the  result  of  PFFM  applied  to  a  pursuit-evasion  game 
in  the  plane  with  control  constraints.  Using  reduced  coordinates  (see  [17]  for 
details)  we  can  re-cast  the  problem  as  an  HJI  equation  by  choosing  in  (3) 

f(x,y,a,b)=vPa-vEb ,  vP  =  2,  vE  =  1 ,  T  =  {\x  -  y\  <  e}  ,  £>0, 

a  G  A  =  L?(0, 1),  b  G  B  =  {(cos #,  sin#)  :  0  G  [7r/4,  77t/4] }  . 
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Figure  7:  Pursuit-Evasion  game  with  control  constraints,  (a)  level  set  of  the 
solution,  (b)  a  couple  of  optimal  trajectories. 

4.3  Patchy  Domain  Decomposition  (PDD)  for  a  class  of 
Hamilton- Jacobi-Bellman  equations 

In  the  spirit  of  the  Domain  Decomposition  technique  described  in  Section  2.3, 
we  developed  a  new  method  for  the  parallelization  of  algorithms  solving  the  HJB 
equation  (1).  The  new  method  is  inspired  by  some  works  of  Ancona  and  Bressan 
(see  [1,  2,  3])  related  to  the  problem  of  the  asymptotic  stabilization  of  a  control 
system.  In  these  works  the  authors  introduced  the  “patchy”  technique,  based  on 
the  division  of  the  domain  in  regions  (patches)  which  are  invariant  with  respect 
to  the  optimal  feedback  dynamics  of  the  problem.  Keeping  this  idea  in  mind,  we 
ended  up  with  a  dynamical  construction  of  a  domain  decomposition  that  exploits 
the  invariance  of  the  sub-domains,  assigning  each  patch  to  a  processor.  The 
resulting  parallel  algorithm  is  very  efficient,  since  it  avoids  the  communication 
costs  between  processors. 

The  main  steps  of  the  construction  are  the  following:  we  first  solve  the  equa¬ 
tion  on  a  coarse  grid  by  means  of  the  classical  (static)  domain  decomposition 
technique.  Then  we  compute  an  approximation  of  the  optimal  feedback  control 
on  the  actual  (fine)  grid  using  in  (2)  a  simple  linear  interpolation  of  the  coarse 
solution.  Finally,  we  define  a  partition  of  the  target  set  and  we  let  flow  each 
part  according  to  the  approximate  optimal  dynamics.  In  practice  each  patch 
is  an  approximation  of  a  tube  of  characteristic  lines  emanating  from  the  target 
set.  We  refer  the  reader  to  [P3]  for  a  detailed  description  of  the  patchy  domain 
decomposition. 

In  Figure  8  we  show  the  invariant  patchy  decompositions  obtained  by  our  algo¬ 
rithm  applied  to  the  following  classical  problems: 

(a)  Eikonal  dynamics  f(x,  y,a)  =  a,  A  =  B( 0, 1)  and  T  =  B( 0,  0.5). 

(b)  Fan  dynamics  f(x,  y,  a)  =  (\x  +  y |  +  0.1)a,  A  =  B( 0, 1)  and  T  =  {x  =  0}. 

(c)  Zermelo  dynamics  f(x,  y ,  a)  =  2.1  a+(2,  0),  A  =  B( 0, 1)  and  T  =  B( 0,  0.5). 
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(d)  Brockett  dynamics  f(x,y,z,a)  =  {a\,a ',2,xa,2  —  ya\),  A  =  [— 5,5]2  and 
T=B(  0,0.25). 


(b) 


n 


(d) 

Figure  8:  Patchy  decompositions  (with  8  patches)  for  some  classical  problems: 
(a)  Eikonal,  (b)  Fan,  (c)  Zermelo,  (d)  Brockett.  Optimal  controls  are  overlaid 
for  the  first  three  problems. 

We  remark  that  the  initial  partition  of  the  target  could  not  be  performed  effi¬ 
ciently  (or  at  all)  for  small  sets,  such  as  lines  or  points.  In  these  cases  we  can 
always  start  our  construction  with  a  single  patch  and  run  the  algorithm  in  serial 
mode  until  the  patch  reaches  a  reasonable  size.  Then  the  solution  is  computed 
in  the  patch  and  it  is  used  as  a  new  target  to  be  divided,  in  order  to  restart  the 
algorithm  in  parallel  mode. 

Similarly,  in  some  problems  the  dynamics  could  induce  very  different  sizes  in 
the  patches,  resulting  in  a  loss  of  performance  for  the  parallel  algorithm  if  one 
processor  finishes  his  job  before  the  others  and  there  are  no  more  patches  in 
the  cue  to  be  assigned.  This  situation  can  be  avoided  by  mixing  iteratively  the 
construction  of  the  patches  with  the  parallel  computation  of  the  solution:  we 
first  run  the  construction  until  the  sizes  of  patches  remain  balanced,  then  we 
compute  in  parallel  the  solution  in  the  patches  and  finally  we  restart  the  proce¬ 
dure  considering  the  computed  solution  as  a  new  target. 

We  also  remark  that  our  method  could  be  applied  in  problems  without  target. 
Indeed,  using  the  coarse  solution,  we  can  approximately  localize  the  sources  of 
characteristics,  where  we  can  put  fictitious  targets  and  then  start  the  patchy 
decomposition. 
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The  performance  of  the  patchy  domain  decomposition  is  very  satisfactory. 
Despite  the  additional  computation  needed  to  build  the  patches,  the  overall  CPU 
time  for  computing  the  solution  is,  for  all  the  test  we  performed,  better  than  that 
of  the  classical  domain  decomposition,  even  better  as  we  increase  the  number  of 
patches.  This  is  due  not  only  to  the  absence  of  transmission  conditions  between 
processors,  but  also  to  a  better  distribution  of  the  computational  resources. 
Indeed,  each  processor  never  remains  idle,  when  it  finishes  his  job  it  can  be 
immediately  re-assigned  to  the  computation  in  another  patch.  Moreover,  if  the 
patches  are  balanced  in  size,  all  the  processors  will  actually  stop  at  the  same 
time.  Nevertheless,  comparing  our  solution  and  the  solution  of  the  classical 
domain  decomposition,  an  additional  error  in  the  approximation  appears,  since 
the  underlying  structured  grid  cannot  accurately  follow  the  characteristic  lines 
defining  the  boundary  of  the  patches.  To  be  precise,  the  “real”  patches  exist  at 
the  continuous  level  and  what  we  compute  is  just  their  projection  on  the  grid. 
Figure  9  shows  the  “patchy  error”  of  the  solutions  for  the  dynamics  (a),(b),(c) 
defined  above. 


Figure  9:  Patchy  error  for  some  classical  problems:  (a)  Eikonal,  (b)  Fan,  (c) 
Zermelo. 


Surprisingly,  the  patchy  error  remains  localized  on  the  boundaries  of  the  patches 
and,  by  numerical  evidence,  it  decreases  as  the  space  step  of  the  grid  goes  to  zero. 
The  theoretical  proof  of  an  error  estimate  of  this  kind  is  still  under  investigation. 


4.4  Features  and  limitations  of  Local  Single  Pass  schemes 
solving  general  Hamilton- Jacobi  equations 

Extending  Fast  Marching  techniques  to  general  Hamilton- Jacobi  equations  is 
still  a  challenging  topic.  In  the  context  of  control  problems  many  numerical 
schemes  were  analyzed  and  employed  (e.g.  finite  difference,  semi-Lagrangian, 
discontinuous  Galerkin  and  finite  volumes)  to  produce  fast  methods  that  can  be 
competitive  with  industrial  solvers  based  on  Pontryagin’s  principle. 

We  recall  that  the  ultimate  goal  is  to  overcome  the  well  known  “curse  of  dimen¬ 
sionality”  related  to  the  Bellman’s  approach  (Dynamic  Programming  Principle) 
for  the  solution  of  such  problems.  In  high  dimension,  terribly  huge  computa¬ 
tional  efforts  are  the  price  one  has  to  pay  in  order  to  easily  get  feedback  controls 
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and  optimal  trajectories,  against  the  open- loop  controls  and  sub-optimal  tra¬ 
jectories  given  by  Pont ryagin- like  methods. 

The  common  philosophy  shared  by  Fast  Marching  methods  is  always  to  reduce 
memory  allocation  and  CPU  times,  by  means  of  local  single-pass  schemes  which 
aim  to 

•  keep  the  computational  efforts  strictly  localized  on  regions  which  are  con¬ 
siderably  small  with  respect  to  the  huge  grids  needed  for  the  numerical 
discretization  of  the  corresponding  equations; 

•  avoiding  useless  computations  minimizing  the  number  of  times  each  grid 
point  is  re-computed. 

Beside  the  development  of  such  fast  solvers,  including  the  methods  described  in 
the  previous  sections,  we  considered  a  deeper  study  [P4]  of  their  features  and 
limitations,  aiming  at  answering  the  following  question: 

can  local  single-pass  methods  solve  every  Hamilton- Jacobi  equation? 

To  this  end  we  reviewed  the  most  known  (and  efficient)  methods  falling  in  the 
class  of  Fast  Marching  solvers,  trying  to  identify  a  comprehensive  classification 
of  the  equations  that  they  can  actually  solve. 

Since  the  introduction  of  the  classical  Fast  Marching  Method  (FMM)  for  the 
Eikonal  equation,  Hamilton- Jacobi  equations  have  been  divided  in  two  classes: 
given  a  grid,  we  have 

(E)  Eikonal-like  equations,  whose  characteristic  lines  coincide  or  lie  in  the 
same  simplex  of  the  gradient  lines  of  their  solutions. 

(!E)  Non  eikonal-like  equations,  for  which  there  exists  at  least  a  point  where 
the  characteristic  line  and  the  gradient  of  the  solution  diverge  enough  to 
not  he  in  the  same  simplex. 


By  construction  (see  Section  2.2)  FMM  clearly  works  for  equations  of  type  (E) 
and  fails  for  equations  of  type  (!E). 

In  the  context  of  Hamilton- Jacobi-Bellman  equations  associated  to  minimum 
time  problems  with  target,  we  ended  up  with  another  possible  classification: 

(D)  Equations  with  smooth  characteristics.  Information  spreads  from  the  tar¬ 
get  to  the  rest  of  the  space  along  smooth  lines,  without  shocks.  The 
corresponding  solution  is  differentiable. 

(!D)  Equations  with  nonsmooth  characteristics.  Information  starts  from  the 
target  and  then  crashes,  creating  shocks.  The  corresponding  solution  is 
Lipschitz  continuous. 
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We  point  out  that  the  original  FMM  was  proposed  for  an  upwind  finite  differ¬ 
ence  discretization  of  the  Eikonal  equation.  In  particular,  a  safe  implementation 
of  the  algorithm  forces  every  node  in  NB  to  be  computed  by  means  of  nodes 
in  ACC  only,  possibly  producing  wrong  values.  Nevertheless,  the  node  in  NB 
with  minimal  value  turns  out  to  be  the  correct  one,  not  affected  by  any  other 
node  in  NB ,  hence  truly  dependent  on  nodes  in  ACC  only. 

On  the  other  hand,  the  semi-Lagrangian  schemes  we  always  employ  for  solving 
Hamilton- Jacobi  equations  (which  are  more  accurate  due  to  their  ability  to  fol¬ 
low  directions  not  aligned  to  the  grid)  usually  compute  nodes  in  NB  by  using 
nodes  in  ACC  and  NB.  In  the  end,  the  minimal  value  rule  among  nodes  in  NB 
still  produces  the  correct  solution  and  this  implies  that  also  in  this  case  the 
correct  node  actually  depends  on  nodes  in  ACC  only. 

This  feature  is  crucial  to  understand  a  point  of  view  that  is  (to  our  best  knowl¬ 
edge)  completely  new  in  the  community  of  Fast  Marching  methods:  in  general, 
the  minimal  value  rule  is  not  an  essential  request  to  make  the  algorithm  work, 
rather  it  could  be  an  obstacle  in  cases  with  strong  anisotropy. 

As  remarked  in  Section  2.2,  accepting  the  node  in  NB  with  minimal  value  cor¬ 
responds  to  compute  the  solution  following  its  level  sets,  a  strategy  that  fails 
for  equations  belonging  to  the  class  (!E)  (e.g.  the  anisotropic  Eikonal  equation). 
Thus  FMM  can  only  solve  equations  in  the  classes  (E  &  D)  or  (E  &  !D). 

On  the  contrary,  equations  falling  in  the  class  (!E  &  D)  (but  clearly  also  equa¬ 
tions  of  type  (E  &  D))  can  be  correctly  solved,  even  faster  than  FMM,  by  simply 
accepting  all  the  nodes  in  NB  that  depend  on  nodes  in  ACC  only.  Indeed,  due 
to  the  assumption  of  being  in  the  class  (D)  (i.e.  there  is  no  shock  and  then  no 
crossing  of  information),  these  nodes  have  values  that  depend  on  information 
coming  from  one  side  only,  driven  by  a  smooth  characteristic.  They  are  com¬ 
puted  starting  from  nodes  in  ACC  belonging  to  the  target  T,  hence  they  are 
correct  by  induction. 

Unfortunately,  the  most  interesting  (but  also  the  most  simple)  equations  as¬ 
sociated  to  control  problems  are  known  to  develop  shocks,  thus  falling  in  the 
class  (!D),  among  which  equations  of  type  (!E  &  !D)  are  even  more  interesting. 
In  these  cases  the  minimal  value  rule  becomes  essential  to  properly  localize  the 
shocks.  Indeed,  the  fact  that  the  narrow  band  NB  expands  as  an  approximation 
of  a  level  set  of  the  solution  implies  that  a  grid  node  close  to  a  shock  is  reached 
from  two  different  directions  at  the  same  time  and  its  value  is  computed  using 
correct  information,  possibly  coming  from  both  sides. 

On  the  other  hand,  if  NB  is  not  an  approximation  of  a  level  set  (a  necessary 
condition,  as  explained  above,  to  deal  with  (!E)  equations),  then  a  grid  node 
close  to  a  shock  could  be  reached  by  NB  at  two  different  times.  In  this  situa¬ 
tion  it  is  impossible  perceiving,  using  local  information  only,  if  the  first  arrival 
is  the  correct  one.  Whenever  an  error  occurs,  it  is  propagated  and  it  cannot 
be  redressed  in  the  future,  since  the  algorithm  continues  the  enlargement  of  NB 
and  ACC  and  never  re-computes  already  accepted  nodes. 

Our  conclusion,  which  is  the  main  achievement  of  this  study,  is  that  in  gen¬ 
eral  equations  of  type  (!E  &  !D)  cannot  be  solved  by  local  single-pass  schemes. 
One  has  to  either  enlarge  the  width  of  the  narrow  band  NB  giving  the  local- 
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ity  up,  and/or  re-compute  the  already  accepted  nodes  giving  the  single-pass 
property  up.  In  both  cases  the  gain  of  Fast  Marching  methods  with  respect  to 
memory  allocation  and  efficiency  is  completely  lost.  In  the  worst  situations,  the 
narrow  band  NB  can  fill  the  whole  domain  and/or  the  number  of  times  one  node 
is  re-computed  can  diverge,  so  that  in  general  the  convergence  of  the  algorithm 
takes  place  only  in  the  limit.  In  the  end  we  recover  the  classical  fixed-point  full 
grid  methods,  that  surely  can  solve  every  equation,  even  if  they  are  terribly  slow 
to  be  applied  to  real  problems  in  high  dimension. 

To  further  convince  the  reader  that  equations  of  type  (!E  &  !D)  are  the  theoret¬ 
ical  limit  for  the  application  of  local  single-pass  methods,  we  finally  remark  that 
we  found  very  pathological  examples  in  our  experiments.  By  coupling  strong 
anisotropy  and  space  inhomogeneity,  we  ended  up  with  two  kinds  of  (!E  &  !D) 
equations  for  which  a  local  single  pass  scheme  cannot  even  start,  not  only  for 
the  problem  of  localizing  the  shocks  described  above. 

In  the  first  case,  there  exists  no  node  in  NB  that  truly  depends  on  nodes  in 
ACC  only:  the  discretization  of  the  equation  introduces,  no  matter  how  fine  is 
the  numerical  grid,  an  artificial  loop-dependency  between  nodes  in  NB  that  can 
be  only  solved  by  means  of  fixed-point  iterations,  thus  breaking  the  single-pass 
property. 

Even  worst,  in  the  second  case,  despite  the  narrow  band  is  stabilized  through 
fixed-point  iterations,  we  still  do  not  find  exact  values  in  NB  (“exact”  compared 
to  solution  of  a  full  grid  method  based  on  the  same  numerical  semi-Lagrangian 
scheme):  the  numerical  grid  introduces  a  kind  of  spurious  reflux  of  information 
from  the  FAR  region  that  has  to  be  taken  into  account  by  enlarging  the  width 
of  the  narrow  band,  thus  breaking  the  locality  property  of  the  method.  Some¬ 
times  this  reflux  can  be  reduced  by  decreasing  the  space  step  of  the  grid,  but  in 
general  it  could  disappear  only  in  the  limit. 

4.5  Numerical  approximation  of  Nash  equilibria 

The  notion  of  Nash  equilibrium  is  considered  one  of  the  most  important  achieve¬ 
ments  in  the  last  century  for  its  impact  on  the  analysis  of  economic  and  social 
phenomena.  While  there  is  an  extensive  literature  dealing  with  Nash  equilibria 
for  static  non-cooperative  games,  the  analysis  of  Nash  equilibria  for  nonzero- 
sum  multi-player  differential  games  is  more  limited  (see  e.g.  the  monographs 
[6,  19]  for  an  presentation  of  this  theory).  Moreover,  only  few  papers  give  a 
characterization  of  the  value  functions  for  the  players  in  terms  of  partial  dif¬ 
ferential  equations,  as  it  was  the  case  for  control  problems  and  for  zero-sum 
differential  games  described  for  example  in  [5].  More  precisely,  we  know  that 
under  suitable  regularity  assumptions,  if  the  value  functions  of  a  non-cooperative 
nonzero-sum  multi-player  differential  game  exist,  they  satisfy  a  system  of  first- 
order  Hamilton- Jacobi  equations,  see  [4].  Typically  theoretical  results  about 
these  problems  are  valid  only  in  very  special  cases  and  they  are  essentially  lim¬ 
ited  to  games  in  one  dimension  and  with  a  simple  dynamics,  see  e.g.  [8,  9,  10]. 
More  important,  it  is  well  known  that  the  system  of  Hamilton- Jacobi  equations 
can  be  ill-posed  and  that  there  are  no  theoretical  results  if  the  dimension  of  the 
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problem  is  greater  than  one. 

From  the  numerical  point  of  view  the  situation  is  even  more  deceiving  since 
few  results  are  available  for  Nash  equilibria  in  the  framework  of  non-cooperative 
games.  In  particular,  we  mention  the  paper  [15]  where  an  approximation  of 
(static)  Nash  equilibria  based  on  Newton  methods  is  presented,  and  the  paper 
[16]  where  the  approximation  is  obtained  via  penalty  methods. 

In  this  context  we  proposed  an  approximation  scheme  for  Nash  equilibria 
of  non-cooperative  differential  games.  Following  Bressan  [8],  we  considered  the 
system  of  stationary  Hamilton- Jacobi  equations  giving  a  characterization  of  the 
value  functions  for  a  class  of  infinite-horizon  games  with  nonlinear  costs  expo¬ 
nentially  discounted  in  time.  We  extended  to  this  system  the  class  of  dynamic 
programming  approximation  schemes  studied  for  zero-sum  differential  games 
(see  [17]  for  a  detailed  analysis  of  the  schemes  corresponding  to  pursuit-evasion 
games  as  well  as  some  numerical  tests).  To  our  best  knowledge,  this  is  the  first 
approximation  scheme  for  computing  Nash  equilibria  in  the  framework  of  differ¬ 
ential  games.  Despite  the  lack  of  theoretical  results,  several  tests  showed  that 
our  scheme  is  convergent  for  some  two-player  games  in  one  and  two  dimensions. 
In  particular,  in  dimension  one,  we  found  that  our  numerical  solution  to  the 
system  of  Hamilton- Jacobi  equations  gives  an  approximation  of  the  “admissible 
solution”  characterized  by  Bressan.  We  refer  the  reader  to  [P5]  for  details  on 
the  new  method. 

4.6  Coupling  Value  Iteration  and  Policy  Iteration  for  solv¬ 
ing  Hamilton- Jacobi  equations 

The  classical  policy  iteration  method  for  control  problems  is  known  to  have  su- 
perlinear  convergence  in  many  interesting  cases,  provided  the  initial  condition  is 
sufficiently  close  to  the  solution  (see  the  monograph  [21]  and  the  precise  result 
in  [23]).  This  limitation  reduces  the  impact  of  the  method  and  forces  to  use 
value  iteration,  which  is  more  expensive  in  terms  of  CPU  time  but  is  globally 
convergent  to  the  optimal  solution  under  rather  general  assumptions. 

We  have  developed  a  new  method  [P6]  for  the  solution  of  Hamilton- Jacobi  equa¬ 
tions  related  to  optimal  control  problems  and  games  that  tries  to  get  the  ad¬ 
vantages  of  both  method  with  an  efficient  coupling.  The  main  idea  is  to  start 
the  computation  using  the  value  iteration  method  and  then  switch  to  the  policy 
iteration  method  when  a  threshold  fixed  on  the  error  is  reached.  The  delicate 
point  is  to  accurately  determine  this  threshold,  in  order  to  avoid  cumbersome 
computation  due  to  the  value  iteration  and,  at  the  same  time,  to  be  reason¬ 
ably  sure  that  the  policy  iteration  method  will  finally  converge  to  the  optimal 
solution.  We  analyzed  several  criteria  for  an  efficient  coupling  in  a  number  of 
examples  in  dimension  two,  three  and  four  illustrating  the  properties  of  the  new 
method. 
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5  Collaborators  and  contacts 

Simone  Cacace  (post-doc  paid  by  the  AFOSR  grant),  Emiliano  Cristiani  (re¬ 
searcher),  Athena  Picarelli  (PhD  student),  Dante  Kalise  (post-doc),  Alessandro 
Alla  (Phd  student). 

5.1  AFRL  Points  of  Contact 

Starting  from  the  2010  review  meeting,  held  in  in  Washington  DC  (August  9-12, 
2010),  we  had  some  contacts  with  Richard  Erwin  (richard.erwin@kirtland.af.mil) 
and  with  Prof.  Meir  Patcher  (AFIT).  R.  Erwin  visited  our  department  on 
September  2011  to  discuss  some  research  topics  relate  to  pursuit-evasion  games. 

6  Workshops  and  conferences 

Beside  the  annual  review  meetings  at  Arlington,  the  results  obtained  during  this 
project  have  been  presented  in  several  invited  talks  at  international  conferences: 

1.  1st  Workshop  on  Computational  Issues  in  Nonlinear  Control,  Monterey, 
November  2009 

2.  EPSRC  Symposium  on  Game  Theory  for  Finance,  Social  and  Biological 
Sciences,  Warwick,  April  2010 

3.  ICIAM  2011  -  Minisymposium  “Efficient  Numerical  Methods  For  Hamilton- 
Jacobi  Equations”,  Organizers:  O.  Bokanowski  and  M.  Falcone,  Vancou¬ 
ver,  July  2011 

4.  SIAM  Conference  on  Control  and  its  Applications,  Baltimore,  July  2011 

5.  IFAC  2011,  Milan,  September  2011 

6.  SADCO  Workshop  “Optimal  Control”,  London,  September,  2011 

7.  Workshop  “Control  and  Optimization  of  PDEs”,  Graz,  October,  2011 

8.  2nd  Workshop  on  Computational  Issues  in  Nonlinear  Control,  Monterey, 
November  2011 

9.  SADCO  Spring  School  &  Workshop  “Applied  and  Numerical  Optimal 
Control”,  Paris,  April  2012 

10.  Minisymposium  at  ORCOS  2012  “Optimal  Control,  Dynamic  Games  and 
Nonlinear  Dynamics”  Vienna  University  of  Technology,  June,  2012. 
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7  Publications 

Here  we  list  the  publication  produced  during  the  three  years  of  the  AFOSR 
project. 

[PI]  M.  Breuss,  E.  Cristiani,  P.  Gwosdek  and  O.  Vogel  An  adaptive 
domain- decomposition  technique  for  parallelization  of  the  Fast  Marching  method , 
Appl.  Math.  Comput.,  218  (2011),  pp.  32-44 

[P2]  S.  Cacace,  E.  Cristiani,  and  M.  Falcone,  A  local  ordered  upwind 
method  for  Hamilton- Jacobi  and  Isaacs  equations ,  Proceedings  of  18th  IFAC 
World  Congress  2011,  pp.  6800-6805. 

[P3]  S.  Cacace,  E.  Cristiani,  M.  Falcone  and  A.  Picarelli,  A  patchy 
dynamic  programmin  scheme  for  a  class  of  Hamilton- Jacobi- Bellman  equations , 
SIAM  J.  Sci.  Comput.,  34  No.  5  (2012),  pp  2625-2649. 

[P4]  S.  Cacace,  E.  Cristiani  and  M.  Falcone,  User  guide  for  single-pass 
methods  solving  Hamilton- Jacobi  equations ,  2012,  in  preparation. 

[P5]  S.  Cacace,  E.  Cristiani  and  M.  Falcone,  Numerical  approximation  of 
Nash  equilibria  for  a  class  of  non-cooperative  differential  games ,  L.  A.  Petrosjan 
e  V.  V.  Mazalov  (Eds.),  Game  Theory  and  Applications,  vol.  16,  Nova  Science 
Publishers,  New  York. 

[P6]  A.  Alla,  M.  Falcone  and  D.  Kalise,  An  accelerated  Policy  Iteration 
algorithm  for  dynamic  programming ,  2012,  in  preparation. 
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